home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / Scientist's Helper src / s.helper.3 / reggauss.c < prev    next >
C/C++ Source or Header  |  1985-11-30  |  11KB  |  384 lines
  1. #include "all.h"
  2. #include "regtabext.h"
  3.  
  4. multifit()
  5. {
  6.     int indList[32], nInd, i, j, k, n, inWord, yCol, zCol, anyNaN, nstore, ierror, counts, flag;
  7.     char c, str[80];
  8.     float t, x[33], y, z, vec[33], *a;
  9.     if( (strcmp(command.cmdWord[1],"type")!=0) &&
  10.          (strcmp(command.cmdWord[1],"keep")!=0) &&
  11.          (strcmp(command.cmdWord[1],"remove")!=0) &&
  12.          (strcmp(command.cmdWord[1],"compute")!=0)  ) {
  13.              ErrMsg(badModifier);
  14.         }
  15.     /*decode independent column list*/
  16.     n = strlen( command.cmdWord[2] );
  17.     command.cmdWord[2][n]=' ';
  18.     command.cmdWord[2][n+1]='\0';
  19.     n++;
  20.     nInd = 0;
  21.     inWord = FALSE;
  22.     for( i=0; i<n; i++ ) {
  23.         c = command.cmdWord[2][i];
  24.         if( (c=='\t') || (c==',') ) c=' ';
  25.         if( (!inWord) && (c!=' ') ) {
  26.             inWord=TRUE;
  27.             nInd++;
  28.             if( nInd>32 ) ErrMsg("too many independent columns");
  29.             j=0;
  30.             str[0]=c;
  31.             }
  32.         else if( inWord &&  (c==' ')  ) {
  33.             j++;
  34.             str[j]='\0';
  35.             SToI( str, &(indList[nInd]) );
  36.             if( GoodCol(indList[nInd])!=0 ) ErrMsg("bad independent column");
  37.             inWord=FALSE;
  38.             }
  39.         else if( inWord && (c!=' ') ) {
  40.             j++;
  41.             str[j]=c;
  42.             } /*end if*/
  43.         } /*end for i*/
  44.     if (nInd<1) ErrMsg("too few independent columns");
  45.         
  46.     SToI( command.cmdWord[3], &yCol );
  47.     if( GoodCol(yCol)!=0 ) ErrMsg("bad dependent column");
  48.     
  49.     /*allocate matrix and zero it*/
  50.     nstore = nInd+1;
  51.     a = NewPtr( 4L * (long)(nstore*nstore) );
  52.     if( MemError()!=noErr ) ErrMsg("not enough memory");
  53.     for( i=0; i<nstore; i++ ) {
  54.         vec[i]=0.0;
  55.         for( j=0; j<nstore; j++ ) {
  56.             *(a+i*nstore+j)=0.0;
  57.             } /*end for j*/
  58.         } /*end for i*/
  59.     
  60.     counts = 0;
  61.     for( i=1; i<=table.header.rows; i++ ) {
  62.         x[0]=1.0;
  63.         anyNaN = FALSE;
  64.         for( j=1; j<=nInd; j++ ) {
  65.             GetTable( i, indList[j], &(x[j]) );
  66.             if( NaN(&x[j]) ) {
  67.                 anyNaN = TRUE;
  68.                 break;
  69.                 } /*end if*/
  70.             } /*end for j*/
  71.         GetTable( i, yCol, &y );
  72.         if( anyNaN || NaN(&y) ) break;
  73.         counts++;
  74.         for( j=0; j<nstore; j++ ) {
  75.             vec[j] += y*x[j];
  76.             for( k=0; k<nstore; k++ ) {
  77.                 *(a+j*nstore+k) += x[j]*x[k];
  78.                 } /*end for k*/
  79.             } /*end for j*/
  80.         } /*end for i*/
  81.             
  82.     gauss(a,vec,nstore,nstore,1.0e-6,&ierror,TRUE); /*solve normal equations*/
  83.     if( ierror!=0 ) {
  84.         WriteLine( "Warning! singular matrix detected" );
  85.         }
  86.     DisposPtr( a );
  87.         
  88.     if (strcmp(command.cmdWord[1],"type")==0) {
  89.         WriteLine("constant term");
  90.         RToS( vec[0], str );
  91.         WriteLine( str );
  92.         WriteLine("coefficients of the indepencent columns");
  93.         for( i=1; i<=nInd; i++ ) {
  94.             IToS( indList[i], str );
  95.             WritePhrase( str );
  96.             WritePhrase( "  " );
  97.             RToS( vec[i], str );
  98.             WriteLine( str );
  99.             } /*end for i*/
  100.         } /*end if*/
  101.     else if( (strcmp(command.cmdWord[1],"keep")==0)  ||
  102.              (strcmp(command.cmdWord[1],"remove")==0)  ) {
  103.         if(strcmp(command.cmdWord[1],"remove")==0) {
  104.             flag=TRUE;
  105.             }
  106.         else {
  107.             flag=FALSE;
  108.             }
  109.         SToI( command.cmdWord[4], &zCol );
  110.         for( i=1; i<=table.header.rows; i++ ) {
  111.             z=vec[0];
  112.             for( j=1; j<=nInd; j++ )  {
  113.                 GetTable( i, indList[j], &t );
  114.                 z+=vec[j]*t;
  115.                 } /*end for j*/
  116.             if (flag) {
  117.                 GetTable( i, yCol, &y );
  118.                 z = y-z;
  119.                 } /*end if*/
  120.             SetTable( i, zCol, z, FALSE );
  121.             } /*end for i*/
  122.         } /*end if*/
  123.         
  124.     /*save coefficients, etc*/
  125.     nCoeffs=nstore;
  126.     for( i=0; i<=nstore; i++ ) coeffs[i]=vec[i];
  127.     IToS( counts, str );
  128.     if( !SetVar( "counts", str ) ) {
  129.         ErrMsg("couldnt set variable");
  130.         }
  131. }
  132.  
  133. noise()
  134. {
  135.     float a, d, r;
  136.     double xnse();
  137.     int xCol, i;
  138.     SToR( command.cmdWord[1], &a, FALSE );
  139.     SToR( command.cmdWord[2], &d, FALSE );
  140.     if( d<=0.0 ) ErrMsg("std dev cant be non-positive");
  141.     SToI( command.cmdWord[3], &xCol );
  142.     for( i=1; i<=table.header.rows; i++ ) {
  143.         r = (float) xnse(a,d);
  144.         SetTable( i, xCol, r, FALSE );
  145.         }
  146. }
  147.  
  148. double xnse(a,d)
  149. float a, d;
  150. /* returns a random number in a sequence with mean a and standard deviation d */
  151. {
  152.     double sum, t, ss;
  153.     int i, n=10;
  154.     sum = 0.0;
  155.     t = sqrt( (double) (12.0*d*d/n) )/65536.0;
  156.     for( i=0; i<n; i++ ) {
  157.         sum +=  ran();
  158.         }
  159.     return( (sum*t + (double)a) );
  160. }
  161.  
  162. polyfit()
  163. {
  164.     int order, i, j, k, xCol, yCol, zCol, counts, ierror, flag;
  165.     float a[7][7], vec[7], xn[14], x, y, z, t;
  166.     char str[40];
  167.     
  168.     SToI( command.cmdWord[1], &order );
  169.     if( (order<1) || (order>6) ) {
  170.         ErrMsg("order must be in range 1-6");
  171.         }
  172.     if( (strcmp(command.cmdWord[2],"type")!=0) &&
  173.          (strcmp(command.cmdWord[2],"keep")!=0) &&
  174.          (strcmp(command.cmdWord[2],"remove")!=0) &&
  175.          (strcmp(command.cmdWord[2],"compute")!=0)  ) {
  176.              ErrMsg(badModifier);
  177.         }
  178.     SToI( command.cmdWord[3], &xCol );
  179.     SToI( command.cmdWord[4], &yCol );
  180.     if( (GoodCol(xCol)!=0) || (GoodCol(yCol)!=0) ) {
  181.         ErrMsg( "bad input column");
  182.         }
  183.         
  184.     for( j=0; j<=order; j++ ) { /*zero matrix and vector*/
  185.             vec[j] = 0.0;
  186.             for( k=0; k<=order; k++ ) {
  187.                 a[j][k] = 0.0;
  188.                 } /*end for k*/
  189.             } /*end for k*/
  190.  
  191.     counts=0;
  192.     i=1;
  193.     while( NextNotNaN( i, xCol, yCol, &i ) ) {
  194.         counts++;
  195.         GetTable( i, xCol, &x );
  196.         GetTable( i, yCol, &y );
  197.         xn[0]=1.0;  /*compute powers of x*/
  198.         for( j=1; j<=(2*order); j++ )  {
  199.             xn[j] = x * xn[j-1];
  200.             } /*end for j*/
  201.         for( j=0; j<=order; j++ ) { /*zero matrix and vector*/
  202.             vec[j] +=  y*xn[j];
  203.             for( k=0; k<=order; k++ ) {
  204.                 a[j][k] += xn[k+j];
  205.                 } /*end for k*/
  206.             } /*end for k*/
  207.         i++;
  208.         } /*end while*/
  209.         
  210.     gauss(a,vec,(order+1),7,1.0e-6,&ierror,TRUE); /*solve normal equations*/
  211.     if( ierror!=0 ) {
  212.         WriteLine( "Warning! singular matrix detected" );
  213.         }
  214.         
  215.     if (strcmp(command.cmdWord[2],"type")==0) {
  216.         WriteLine("coefficients of x to the n, starting with n=0");
  217.         for( i=0; i<=order; i++ ) {
  218.             RToS( vec[i], str );
  219.             WriteLine( str );
  220.             } /*end for i*/
  221.         } /*end if*/
  222.     else if( (strcmp(command.cmdWord[2],"keep")==0)  ||
  223.              (strcmp(command.cmdWord[2],"remove")==0)  ) {
  224.         if(strcmp(command.cmdWord[2],"remove")==0) {
  225.             flag=TRUE;
  226.             }
  227.         else {
  228.             flag=FALSE;
  229.             }
  230.         SToI( command.cmdWord[5], &zCol );
  231.         for( i=1; i<=table.header.rows; i++ ) {
  232.             GetTable( i, xCol, &x );
  233.             xn[0]=1.0;
  234.             z=vec[0];
  235.             for( j=1; j<=order; j++ )  {
  236.                 xn[j]=x*xn[j-1];
  237.                 z+=xn[j]*vec[j];
  238.                 } /*end for j*/
  239.             if (flag) {
  240.                 GetTable( i, yCol, &y );
  241.                 z = y-z;
  242.                 } /*end if*/
  243.             SetTable( i, zCol, z, FALSE );
  244.             } /*end for i*/
  245.         } /*end if*/
  246.         
  247.     /*save coefficients, etc*/
  248.     nCoeffs=order;
  249.     for( i=0; i<=order; i++ ) coeffs[i]=vec[i];
  250.     IToS( counts, str );
  251.     if( !SetVar( "counts", str ) ) {
  252.         ErrMsg("couldnt set variable");
  253.         }
  254.  
  255. }
  256.  
  257. gauss(a,vec,n,nstore,test,ierror,itriag)
  258. float *a, vec[], test;
  259. int n, nstore, *ierror, itriag;
  260. {
  261.  
  262. /*subroutine gauss, by william menke, july 1978 (modified feb 1983, nov 85) */
  263.  
  264. /* a subroutine to solve a system of n linear equations in n unknowns, where n doesn't exceed 32 */
  265. /*gaussian reduction with partial pivoting is used                 */
  266. /*     a        (sent, destroyed)    n by n matrix                    */
  267. /*    vec        (sent, overwritten)    n vector, replaced with solution    */
  268. /*     nstore    (sent)            dimension of a                    */
  269. /*    test        (sent)            small pos no for div by zero check    */
  270. /*    ierror    (returned)        error code. zero on no error        */
  271. /*    itriag    (sent)            matrix triangularized only on TRUE    */
  272. /*                            useful when solving multiple        */
  273. /*                            systems with same a             */
  274.     static int isub[32], l1;
  275.     int line[32], iet, ieb, i, j, k, l, j2;
  276.     float big, testa, b, sum;
  277.     
  278.  
  279.     iet=0;    /* initial error flags, one for triagularization*/
  280.     ieb=0;    /* one for backsolving */
  281.  
  282. /* triangularize the matrix a, replacing the zero elements of the triangularized matrix */
  283. /* with the coefficients needed to transform the vector vec */
  284.  
  285.     if (itriag) {    /* triangularize matrix */
  286.  
  287.          for( j=0; j<n; j++ ) {      /*line is an array of flags*/
  288.             line[j]=0;   /* elements of a are not moved during pivoting. line=0 flags unused lines */
  289.             }    /*end for j*/
  290.             
  291.         for( j=0; j<n-1; j++ ) { /*  triangularize matrix by partial pivoting */
  292.  
  293.                        big = 0.0;     /* find biggest element in j-th column of unused portion of matrix*/
  294.                for( l1=0; l1<n; l1++ ) {
  295.                                if( line[l1]==0 ) {
  296.                                        testa=(float) fabs( (double) (*(a+l1*nstore+j)) );
  297.                                        if (testa>big) {
  298.                         i=l1;
  299.                         big=testa;
  300.                         } /*end if*/
  301.                     } /*end if*/
  302.                 } /*end for l1*/
  303.                        if( big<=test) {   /* test for div by 0 */
  304.                                iet=1;
  305.                                } /*end if*/
  306.  
  307.                        line[i]=1;  /* selected unused line becomes used line */
  308.                        isub[j]=i;  /* isub points to j-th row of triang. matrix */
  309.  
  310.                        sum=1.0/(*(a+i*nstore+j));  /* reduce matrix towards triangle */
  311.                for( k=0; k<n; k++ ) {
  312.                 if( line[k]==0 ) {
  313.                     b=(*(a+k*nstore+j))*sum;
  314.                            for( l=j+1; l<n; l++ ) {
  315.                                                *(a+k*nstore+l)=(*(a+k*nstore+l))-b*(*(a+i*nstore+l));
  316.                            } /*end for l*/
  317.                                        *(a+k*nstore+j)=b;
  318.                     } /*end if*/
  319.                 } /*end for k*/
  320.  
  321.             } /*end for j*/
  322.  
  323.                for( j=0; j<n; j++ ) { /*find last unused row and set its pointer*/
  324.             /*  this row contians the apex of the triangle*/
  325.             if( line[j]==0) {
  326.                 l1=j;   /*apex of triangle*/
  327.                 isub[n-1]=j;
  328.                 break;
  329.                 } /*end if*/
  330.             } /*end for j*/
  331.  
  332.         } /*end if itriag true*/
  333.         
  334.     /*start backsolving*/
  335.     
  336.     for( i=0; i<n; i++ ) {        /* invert pointers. line(i) now gives row no in triang matrix */
  337.                         /* of i-th row of actual matrix */
  338.         line[isub[i]] = i;
  339.         } /*end for i*/
  340.  
  341.       for( j=0; j<n-1; j++) {  /* transform the vector to match triang. matrix*/
  342.                b=vec[isub[j]];
  343.                for( k=0; k<n; k++ ) {
  344.                       if (line[k]>j) {    /* skip elements outside of triangle*/
  345.                                 vec[k]=vec[k]-(*(a+k*nstore+j))*b;
  346.                 } /*end if*/
  347.             } /*end for k*/
  348.         } /*end for j*/
  349.  
  350.       b = *(a+l1*nstore+(n-1));   /*apex of triangle*/
  351.       if( ((float)fabs( (double) b))<=test) {  /*check for div by zero in backsolving*/
  352.         ieb=2;
  353.         } /*end if*/
  354.       vec[isub[n-1]]=vec[isub[n-1]]/b;
  355.  
  356.       for( j=n-2; j>=0; j-- ) {    /* backsolve rest of triangle*/
  357.         sum=vec[isub[j]];
  358.         for( j2=j+1; j2<n; j2++ ) {
  359.             sum = sum - vec[isub[j2]] * (*(a+isub[j]*nstore+j2));
  360.             } /*end for j2*/
  361.             b = *(a+isub[j]*nstore+j);
  362.                if( ((float)fabs((double)b))<=test) {    /* test for div by 0 in backsolving */
  363.             ieb=2;
  364.             } /*end if*/
  365.         vec[isub[j]]=sum/b;   /*solution returned in vec*/
  366.         } /*end for j*/
  367.  
  368. /*put the solution vector into the proper order*/
  369.  
  370.       for( i=0; i<n; i++ ) {    /* reorder solution */
  371.         for( k=i; k<n; k++ ) {  /* search for i-th solution element */
  372.             if( line[k]==i ) {
  373.                 j=k;
  374.                 break;
  375.                 } /*end if*/
  376.             } /*end for k*/
  377.                b = vec[j];       /* swap solution and pointer elements*/
  378.                vec[j] = vec[i];
  379.                vec[i] = b;
  380.                line[j] = line[i];
  381.         } /*end for i*/
  382.  
  383.       *ierror = iet + ieb;   /* set final error flag*/
  384. }